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Abstract. Colliding hypersonic flows play a decisive role in many astrophysical objects. They contribute, for 
example, to the molecular cloud structure, the X-ray emission of O-stars, differentiation of galactic sheets, ap- 
pearance of wind-driven structures, or, possibly, to the prompt emission of 7-ray bursts. Our intention is thorough 
investigation of the turbulent interaction zone of such flows, the cold dense layer (CDL). In this paper, we focus 
on the idealized model of a 2D plane parallel isothermal slab and on symmetric settings, where both flows have 
equal parameters. We performed a set of high-resolution simulations with upwind Mach-numbers, 5 < M u < 90. 
We find that the CDL is irregularly shaped and has a patchy and filamentary interior. The size of these structures 
increases with £ c di, the extension of the CDL. On average, but not at each moment, the solution is nearly self-similar 
and only depends on M u . We give the corresponding analytical expressions, with numerical constants derived from 
the simulation results. In particular, we find the root-mean-square Mach-number to scale as M rms ~ 0.2A/ U . The 
mean density, p m ~ 30p u is independent of M u . The fraction / c ff of the upwind kinetic energy that survives 
shock passage scales as / e g = 1 — M^ s 6 . This dependence persists if the upwind flow parameters differ from 
one side to the other of the CDL, indicating that the turbulence within the CDL and its driving are mutually 
coupled. Another finding points in the same direction, namely that the auto-correlation length of the confining 
shocks and the characteristic length scale of the turbulence within the CDL are proportional. Larger upstream 
Mach-numbers lead to a faster expanding CDL, confining interfaces that are less inclined with respect to the 
upstream flow direction, more efficient driving, and finer interior structure with respect to the extension of the 
CDL. 
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1. Introduction 

Supersonically turbulent, shock-bound interaction zones 
are important for a variety of astrophysical objects. They 
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and Herbig-Haro objects |]Vlatzner fc McKedll999|) . 

So far, the shape and turbulent interior of shock-bound 
interaction zones have been mostly studied separately. In 
this paper we focus on the system as a whole, stressing 
that upwind flows, confining interfaces of the interaction 
zone, and the interior structure of this zone form a tightly 
coupled system. The turbulence within the interaction 
zone affects the shape of the confining shocks, which in 
turn determines how much energy is thermalized at these 
shocks and how much energy remains available for driving 
the turbulence. 

A variety of papers have been written on the shape 
and stability of 2D interaction zones, of which we men- 
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tion only a few. IVishniacI (^394) shows by analytical 
means that geometrically thin, isothermal, 2D, planar, 
shock-bounded slabs are non-linearly unstable, coining the 
term non- linear thin shell instability, or NTSI, for this 
instability. iBlondin fc Marks! l)l996j) essentially reproduce 
these analytical predictions numerically, also mentioning 
the occurrence of supersonic turbulence within the slab. 
Performing 2D radiative and isotherm a l simu lations of col- 
liding molecular clouds, iKlein et all (^f)98) observe the 
complex shaping and instability of the collision zone. The 
role of a radiat ive cooling layer has been a ddressed by 
several authors. IStrickland fc Blondinl l)l995|) numerically 
investigated flows against a wall in 2D, finding that an 
unstable cooling layer introduces disturbances in the in- 
terface separating the cooling layer from the cooled mat- 
ter. Looking at colliding flow s instead of a flow against a 
wall, Walder & Folini l)l998[) show that one unstable cool- 
ing layer is sufficient to destabilize both confining inter- 
faces of the cooled matter. In addition, the cooled mat- 
ter becomes supersonically turbulent. If self-gravity is in- 
cluded fragmentation of the interaction zone is observed 
l)Anninos fe Normanlll99ei Irlunter et aTlll986|) . 

An overwhelming amount of literature meanwhile 
exists on supersonic turbulence. At least part of 
this attention arises because it is thought that su- 
personic turbulence can explain the structuring and 
support of molecular clouds and thus that it plays 
a decisive role in star formation. A comprehensive 
view of this issue can be found in the recent re- 
views by iMac Low fc Kle ssenl ll2QQ4|). [Elmegreen fc Scald 
(|2004l) . and lScalo fc Elmegreen! ^200i) . Of particular in- 
terest for the work we present here is the paper by 
IMac Low! lll999fl , where Fig. 4 shows that the wave length 
of the driving is apparent in the spatial scale of the 
turbulent structure for monochromatically driven turbu- 
lence in a 3D periodic box. The possible importance of 
the finite size of the slab was recently pointed out by 
iBurkert fc Hartmannl l)2004[h 

We are trying to make four points with this pa- 
per. First, we argue that, within the frame of isother- 
mal Euler equations and in infinite space, the solution 
may be self-similar and dependent only on the upstream 
Mach-number, at least to first approximation. Based on 
this assumption, we give expressions for average quanti- 
ties of the slab. Second, we show that the numerical so- 
lution, which is defined only on a finite computational 
domain and includes (implicit) numerical dissipation, re- 
mains close to self-similar, as long as the width of the slab 
is small and the root-mean-square Mach-number larger 
than one. Third, we stress the tight mutual coupling be- 
tween the turbulence and its driving. Fourth, we point 
out that spatial scales generally grow with extension f c di 
of the interaction zone, but decrease with increasing up- 
stream Mach-number M u . 

Results are based on a set of simulations that differ 
only in their upwind Mach- numbers. In this paper we 
restrict the analysis of these simulations to the above- 
mentioned three objectives. We postpone a more detailed 



analysis of the interior structure of the interaction zone to 
a subsequent paper. 

In the following, we first give the details of our phys- 
ical model and numerical method in Sect. . In Sect. |21 
we derive the self-similar scaling relations. The numerical 
results are present in Sect. 21 Discussion follows in Sect. [SI 
and conclusions in Sect. |H1 

2. Physical model and numerical method 

The numerical treatment of supersonic turbulence is an 
issue in its own right, so we start this section with a brief 
summary of some results that are relevant to the present 
work. We then specify the physical model we consider, 
explain the numerical method we use and the simulations 
we perform. 

2.1. Simulating supersonic turbulence 

The shock-compressed layer studied in this paper is super- 
sonically turbulent with root-mean-square Mach-numbers 
between about 1 and 10. An important fraction of the ki- 
netic energy is dissipated in shocks. Euler equations are 
sufficient for describing this part of the problem. A cas- 
cade transfers the remaining energy to higher and higher 
wave numbers until it is finally destroyed on the viscous 
dissipation scale. To also capture this part of the problem, 
the compressible Navier-Stokes equations should be used; 
however, the range of spatial scales associated with the 
energy cascade exceeds the capacity of any computer by 
far. 

In subsonic turbulence, one way out is to use a suitable 
sub-grid scale model. The model is used to compute an 
effective viscosity coefficient, which should mimic the cas- 
cading between the smallest scale still resolved by the nu- 
merical grid and the viscous dissipation scale as precisely 
as possible. This coefficient is then used in the Navier- 
Stoke s equations instead of the physical viscosity l|Lesieurl 
Il999l) . For the approach to work it is essential that the 
effective viscosity obtained from the sub-grid scale model 
exceeds the (implicit) numerical viscosity of the overall 
numerical scheme. This can be achieved in s ubsonic tur - 
bulence by the use of low-dissipation schemes (|Lelel|l9 92) . 

In supersonic turbulence, explicit sub-grid scale mod- 
eling so far does not exist in the above sense. The ba- 
sic reason is that the numerical treatment of supersonic 
turbulence requires schemes that can treat shocks appro- 
priately, such as the widely used shock capturing schemes. 
The (implicit) numerical viscosity of such schemes is, how- 
ever, much too large to match the abov e requirement, 
even if the schemes are of a high order ( Gamier et alJ 
Il999t IPorter et ai1ll992f) . One strategy for this case, the 
so called MILES approach (mon otone integrated l arge- 
eddy simulation), w as proposed bvlBoris et all l)l992|) and 
further explored by IPorter et all l)l992l \l99^ . The basic 
claim is that the numerical viscosity inherent to shock cap- 
turing schemes llrlirschlll 99,4 iLeVemifJbnolj) acts already 
as a physically correct sub-grid scale model. Solving the 
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Eulcr equations by means of a shock capturing scheme 
thus should yield the correct physical answer. 

The validity of the claim that implicit numerical vis- 
cosity alone leads to a correct ph ysical solution was in- 
vestigated bv lGarnier et al.1 l|l999l) for a selection of shock 
capturing schemes, among them a MUSCL-scheme (mono- 
tone upwind scheme for conservation laws) similar to the 
one we use (see Sect. I2.3)l . For the cases considered (es- 
sentially decaying subsonic), they find that the scheme 
indeed acts as a (very dissipative) sub-grid scale model 
in that it preserves the flow from energy accumulation on 
small spatial scales. However, they also find that struc- 
tures defined on less than 5 gr id points are affect ed by 
substantial numerical damping. IPorter et all il994\ find, 
in addition, that the dissipation properties of their scheme 
(MUSCL with PPM) are highly non-linear, and also they 
depend not only on the grid spacing but also on the wave 
length of the flow structure. Structures on less than 32 
grid points are affected by numerical damping. 

We rely on the MILES approach in this paper for the 
lack of a better model, although, to our knowledge, the va- 
lidity and quality of the approach has never been tested for 
supersonic turbulence. The numerical solutions we obtain 
are thus rather solutions of the Navier-Stokes equations. 
Nevertheless, as dissipation in shocks by far dominates nu- 
merical dissipation, we expect the 'Euler character' of the 
solution to prevail. 
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Fig. 1. Sketch of physical model problem, p^ M;, and si 
denote the density, Mach-number, and confining shock of 
the left (i — I) and right (i = r) flow, p and M denote 
the density and Mach-number of the CDL. a is the abso- 
lute value of the angle between the x-axis and the tangent 
to the shock. CDL is the shock-compressed interaction 
zone. The dashed rectangle indicates the computational 
domain with y-extension Y. Periodic boundary conditions 
in y-direction imply periodic continuation of the solution 
(dotted continuation of left and right shock). 



2.2. The model problem 

The model problem we consider consists of a 2D, plane- 
parallel, infinitely extended, isothermal, shock compressed 
slab. A sketch is given in Fig. ^ Two high Mach-number 
flows, oriented parallel (left flow, subscript I) and anti- 
parallel (right flow, subscript r) to the x-direction, col- 
lide head on. The resulting high-density interaction zone, 
the shock compressed slab, is oriented in the y-direction. 
We denote this interaction zone by CDL for 'cold dense 
la yer' to remain consistent wi th notation used already 
in I Walder fc Folinil I 



1996, 1998). We investigated this sys- 
tem within the frame of Euler equations (but see also 
Sect. 12. l|) . together with a polytropic equation of state, 



dp 
dt 



+ V \ pv <X> v 
at \ p 



P -I 



dE 
~dt 



■ V(« (E + p)) 



0, 


(1) 


0, 


(2) 


0, 


(3) 


p/(7-l). 


(4) 



Here, p is the particle density, p the average mass per par- 
ticle, v — (v x ,v y ) is the velocity vector, p thermal pres- 
sure, / the identity tensor, e the thermal energy density, 
and E = pv 2 /2 + e the total energy density. For the poly- 
tropic exponent, we choose 7 = 1.000001. This value guar- 
antees that jump conditions and wave speeds of a Mach-90 
shock are within 0.01 per cent of the isothermal values. 



Within the frame of this paper we consider only sym- 
metric settings, where the left (subscript I) and right (sub- 
script r) colliding flow have identical parameters (sub- 
script u for upstream): p\ = p r = p u and \v\ \ = \v T \ = v u . 

We look at the problem in a dimensionlcss form and ex- 
press velocities in units of the isothermal sound speed a = 
^/X7cb7m> with T the temperature and fee the Boltzmann 
constant. Densities we express in terms of the upstream 
density p u . Finally, we express lengths in units of Yo, the 
smallest y-extent of the computational domain we used. 
This artificial choice is necessary as there is no natural 
time-independent length scale to the problem (see Sect.|3J). 



2.3. Numerical method 

Our results were with the AMRC ART-code 1 . We 
used the multidimensional high -resolution fin ite-volume- 
integration scheme developed bv lColellal <)l990|) on the ba- 
sis of a Cartesian mesh. Tests showed that this algorithm, 
compared to dimensional splitting schemes, is significantly 
more accurate in capturing flow features not aligned with 



1 AM RC ART is part o f the A-MAZE code- 
package iWalder fc Folinill2000al) . which contains 3D adaptive 
MHD and radiative transfer codes. The package, along with a 
brief description, is publicly available at 
http://www.astro.phys.ethz.ch/staff/folini/folini.html or 
http://www.astro.phys.ethz.ch/staff/walder/walder.html 



4 



Folirii & Walder: Compressible turbulence in shock-bounded interaction zones 



the axis of the mesh. In all our simulations we used a ver- 
sion of the scheme that is (formally) second order accurate 
in space and in time for smooth flows. 

We combine thi s integration sc heme with the adaptive 
mesh algorithm bv iBereerl l)l985|) . While a rather coarse 
mesh was sufficient for the upwind flows, the turbulent 
CDL was resolved on a much finer scale. 

We found it useful to have our CDL moving in positive 
x-direction at a speed of about Mach 20-40. If the CDL 
was essentially stationary with respect to the computa- 
tional grid, we observed alignment effects of strong shocks 
that were nearly parallel to a cell interface (in y-direction). 
Through the global motion of the CDL, which implied su- 
personic motion of the confining shocks with respect to the 
computational grid, we got rid of this problem. We checked 
that this procedure introduced no systematic effects into 
the solution. The problem of alignment effects when 
dealing with high Mach-number flows, nearly stationary 
shocks, and high order upwin d schemes is well known and 
not particular to our scheme llColella fc Woodwardl ll984: 
lOuirkl Il994t I.Tasak fc Welier! Il995h . Other work arounds 
exist, such as smoothing of interfaces by additional vis- 
cosity, which is often applied in PPM implementations. 

2.3.1. Numerical settings and integration time 

In the x-direction, our computational domain extended 
over 200 Yo. The y-extent Y of our domain varied between 
simulations, Yo < Y < 6Yo (see Table |B~T|) . Boundary 
conditions at the left and right boundaries (x-direction) 
were 'supersonic inflow'. In the y-direction we had periodic 
boundary conditions. The cell size at the coarsest level was 
0.2 Yo. The cells at the finest level, covering the CDL, were 
smaller by a factor 2 6 to 2 9 , yielding between 320 and 2560 
cells over a distance Yo (depending on the simulation, see 
Table [ED}. 

As will be shown, the relevant time-dependent quan- 
tity for the evolution of CDL mean quantities is the 
average x-extension of the CDL, £ c di- We defined it as 
4di = V/Y, where V is the 2D volume of the CDL. For 
later use we also introduce the volume integrated den- 
sity 77i c di = / v p, the mean density p m = m c di/V, and 
the average column density N = m c di/Y = /O m ^ c di- The 
last quantity was made dimensionless by division through 
N = p u Yo. We stopped most simulations at £ c di = Y/2. 

2.3.2. Initial conditions 

We investigated three different initial conditions, 1=0,1,2. 

1=0: No CDL exists at t = 0. The left and right flows 
are initially separated by a single interface. The interface 
is wiggled with a single, sinusoidal mode of wave length 
0.1 Y and amplitude 0.0195 Yo (about 3 to 25 grid cells, 
depending on the discretization). 

1=1: A CDL is present at time t = 0. It has a col- 
umn density of N = 14 No and a thickness of 0.03125 Yo- 
The confining shocks are both wiggled, with the same si- 



nusoidal mode and amplitude as the interface in the case 
1=0. The mass within the CDL is at rest and of constant 
density, p — p u M^, the density the CDL would have in 
ID. Note that this initialization implies some violation of 
the Rankine-Hugoniot jump conditions at the interfaces. 

1=2: A CDL is present at time t = 0, with column 
density N = 56 No and a thickness of 0.125 Yo. The right 
shock is wiggled as for 1=1, the left shock is straight. The 
density and velocity in the CDL are set as for 1=1. 

We stress that the initial wiggling of the shocks is not 
compelling. The only effect of this wiggling is to speed up 
the initial phase of the evolution. Test cases using another 
wiggling or starting from straight shocks end up like the 
simulations we are going to present in the following. 

We would like to add a side note on this last point, from 
our observation that the slab is also destabilized when 
bou nd by straight shocks. T his has already been reported 
by iBlondin fc Marks! l|l996h . who ascribe d the destabi- 
lizatio n to 'numerical noise'. Meanwhile, iRobinet et alJ 
(2000) have investigated what is called the carbuncle phe- 
nomenon in some more detail. They showed that - con- 
trary to what has been believed so far - a single straight 
shock is linearly unstable for exactly one mode associated 
to the upstream Mach-number of Af cr i t = [(5 + 7)/(3 — 
7)] 1 / 2 . For isothermal conditions, this yields M cr ;t = -\/3. 
They also showed that this single unstable mode is suffi- 
cient for making straight shocks aligned with the mesh 
numerically unstable at all Mach-numbers if the com- 
putation is done with a low-viscosity, high-order, shock- 
capturing scheme. To what degree this instability for 
a straight shock of any Mach-number is really physical 
seems an open question to us. 

2.4. The different runs 

The runs we performed differ in their upwind Mach- 
numbers, which lie in a range 5 ^ M u <, 90, as well as 
in their initialization, numerical discretization, and the y- 
extent of the domain. The labels of the different runs are 
built up as MJ.R.Y. Here, M is the upwind Mach-number, 
I the initialization (0, 1, or 2), R gives the refinement of 
the spatial discretization, relative to the coarsest grid sim- 
ulation we performed (1, 2, 4, or 8). R=l corresponds to 
a finest cell size of about 3 • 10~ 3 Yo, R=2 indicates a 
twice smaller cell size. Y is the domain size (1, 2, 4, or 
6) in units of Yo- For example, R22_0.2.4 denotes a run 
with M u = 22, initialization 1=0, finest cell size about 
1.5 • 10" 3 Y , and y-extent 4Y . 

The runs we performed are listed in Table IB. II 
Individual columns in Tabic |B~T1 contain (column number 
in square brackets): label of run [1], following the scheme 
label= M U _I.R.Y, where I is the initial condition, R the 
refinement factor such that cell size = 3.125 • 10~ 3 Yo/R, 
and Y is the y-extension of the computational domain in 
units of Yo; Mach-number of upstream flow, M u [2]; stop- 
ping time of simulation in terms of £(N) [3]; y-averaged 
x-extension of CDL at stopping time, relative to y-extent 
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Fig. 2. The self-similar ID solution of isothermal colliding 
supersonic flows in density (top) and velocity (bottom). 
The interaction zone (labeled CDL) is bounded by two 
shocks, si and s r , having speeds vf and vf in the rest 
frame of the CDL. The density and velocity of the ID 
interaction zone, we denote by pid and vid, respectively. 

of computational domain, 4di/Y [4]; average quantities 
[5-9] of: rms Mach-number, M rms [5]; mean density in 
units of upstream density p m /p u [6]; shock length in 
units of y-domain, £ s h/Y [7]; driving efficiency, f e g [8]; 
averages taken over 10 < £(N) < 70 for 1=0 and over 
60 < £(N) < 120 for 1=1, for 1=2 we give the values at 
the end of the simulation in parentheses instead. 

3. Scaling properties of the model problem 

Within the frame of Euler equations and in infinite space, 
the problem of isothermal supersonically colliding flows 
can be solved analytically in ID. The solution, sketched in 
Fig.|21and Sect. 13. ll is self-similar and depends only on two 
free parameters, the Mach-numbers of the left and right 
upwind flow. In 2D th e situation is more complicated: 
the solution is unstable l|Vishniadll99llBlondin fc Marksl 
the shocks confining the CDL are non-stationary 
and oblique, the interior of the CDL is supersonically tur- 
bulent. 

Nevertheless, in infinite space it seems reasonable to 
assume that the solution, on average, may still evolve in a 
self-similar manner. We base this assumption on the fol- 
lowing two observations. First, the isothermal Euler equa- 
tions are scale-free in infinite space. Second, the free pa- 
rameters of the problem (p u , M u , and a) do not introduce 
any fixed length or time scale. Under these conditions, it is 
possible that the solution also does not depend on length 
or time separately, but only on their ratio. If so, all length 
scales should evolve equally with time, which implies, in 



particular, that the solution then should not depend on 
the extension of the CDL. We stress, however, that we 
have no proof of the above assumption of self- similarity. 

In the remainder of this section, we elaborate a bit 
further on the implications of the assumed self-similarity. 
In Sect. 01 we will see that the relations derived here give 
a good approximation of the numerical results, but we 
stress already here three important points. The numerical 
simulations are carried out in finite space (not infinite); 
numerical dissipation might play a role; and the simu- 
lations are stopped for the most part while the CDL is 
still small, about half the size of the y-extent of the com- 
putational domain. Important aspects that can only be 
obtained from the numerical solution include quantities 
related to the driving of the turbulence, the values of pro- 
portionality constants, and the interior structure of the 
CDL. We neglect this last aspect, however, in the current 
paper to focus on mean quantities instead. 

3.1. Self-similar ID solution 

Denoting the density and velocity of the CDL by pid and 
vid , and those of the left and right upwind flows by pi and 
Vi (i = l,r), the solution in the rest frame of the CDL is 
given by 



Pid/ Pi 



Mf 
0, 



1 w Mf 



\v?\ = aMi/{M? - 1) w a/Mi « a. 



(5) 
(6) 
(7) 



Here, if is the velocity of the confining shocks and a is 
again the isothermal sound speed. The approximations 
hold for large Mach-numbers. The self-similar character 
is apparent: the solution is not a function of x and t but 
only a function of x/t through vf . 

A relation between characteristic length and time 
scales of the solution, the self-similarity variable Kid, can 
be obtained as follows. As a length scale, we take the spa- 
tial extension i\d of the CDL, and as a time scale the time 
r needed to accumulate the corresponding column density 
Aid- From the relations 



Aid = PicAd- 
and 

Nld = t (pivi + p r v r ) 

and using p\j ' p T = M^ /My (see Eq. |SJ, we obtain 



lid 

T 



Ml + M r 
1 Ml ■ M r 



(8) 



(9) 



(10) 



Thus for strong shocks Kid is nothing else than |t;f 



Specializing to symmetric settings (1 
Ml and «j ld = 2a/M u . 



r) yields pi d /p u 



3.2. Scaling properties of the 2D symmetric solution 

In the following, we derive scaling relations for the 2D 
solution, assuming self-similarity. We confront these rela- 
tions with corresponding numerical results in Sect.QJ 
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3.2.1. Density, Mach-number, self-similarity variable 

In the following, all velocities are again given in the rest 
frame of the CDL and we assume that a self-similar so- 
lution exists. A natural choice for the (constant) self- 
similarity variable then is again K2d = £ c di/T. Using the 
definitions of Sect. !2.3TTl we must have, as in the ID case, 



N = 2rp u v u . 



(11) 
(12) 



Dividing the two equations through each other yields 
K-2d = 2p u v u /p m . As K2d is a constant, the CDL mean den- 
sity p m must be constant in time. The root-mean-square 
velocity i> 2 ms then has to be constant in time as well, at 
least if the CDL density and velocity, p and u, are uncor- 
rected (in which case we can replace the average over the 
product pv 2 by the product of the averages of p and v 2 ) 
and if kinetic pressure dominates over thermal pressure. 
This can be seen from equating the total upwind pressure 
with the total pressure within the CDL, 



( 2 . 2\ / 2 , 2 \ 

Pu{a +v u )= p m (a +V Ims ). 



(13) 



The simplest ansatz for p m and v Tms is that they only 
depend on the upstream Mach-number, 



p m /Pu = ViM?\ 
Vrms/a = r/ 2 M^ 2 . 



(14) 
(15) 



Using the ansatz for p m we obtain a first expression for 
K2d from Eas. Illland ll2l 

K 2d = 2a V ^M^^ oc aMl~^. (16) 
A second expression for «2d, we obtain from Eg. 1131 

a 2 N, 



p u a (1 + M u ) = p m (a + w rms ) 



cell 



-(1+ V 2 M 2 ^). (17) 



Again using Eq. ^|to replace N, one obtains 



n 2 d = 2aM u 



I + VIM 2 ^ 
1 + Ml 



2ar 1 2 Mf 2 - 1 cx aMf" 1 . (18) 



The approximation is good for high Mach-number flows, 
with rfeMff 2 >> 1, and for /3 2 > 0, which is, however, to 
be expected for supersonic turbulence. Comparing Eas. 1161 
and ^| gives 



Vi 



fa = 1 - ft/2, 
l 



(19) 
(20) 



3.2.2. Driving energy 



From energy conservation, we have fdiss = £drv ~ £kin- 
Here £ d rv is the energy flux density entering the CDL per 
time and per unit length in the y-direction, and £diss de- 
notes the energy density dissipated per time within an 
average column of length £ c d\ of the CDL. Finally, £kin is 
the change per time of the kinetic energy contained within 
such an average column. We first turn to the driving en- 
ergy £drv and come back to fdiss and £kin in Sect. 13.231 



Part of the total (left plus right) upwind kinetic energy 
flux density, .Fe kJn ,u = Pui^J, is thermalized at the shocks 
confining the CDL. The remaining part, £drv, drives the 
turbulence in the CDL. We assume that fdrv and !F eunjU 
are related by a function of the upwind Mach-number only, 



£drv = /effC-MuJ-^ekin.u- 



(21) 



We call the function / e ff the driving efficiency. An expres- 
sion for f e g can be derived by using the jump conditions 
for strong, oblique shocks, 



Pd 
v±,d 



p u Ml u = p u M 2 sin 2 a, 



-2 = a 

M u sina' 
aM u cos a. 



(22) 



The subscript d denotes downstream quantities, right after 
shock passage; the subscripts _L and || denote flow com- 
ponents perpendicular and parallel to the shock, respec- 
tively; and a is given in Fig.Q] Using Eq. [23 we obtain 



6 _ I 

^drv y 



as — - — v± d 



Puvl 
2Y 



1 



dy(l - sin 2 a + ., ,. 

Yi, r Af4sin 2 a 



(23) 



where the integral over si, r and Yi jr runs over both shocks 
and where it was used that sin a ds = dy. The last term 
on the right hand side of Eq. [33] is omitted in the fol- 
lowing. This is justified, as the shocks we observe in our 
simulations fulfill sincv >> M~ 2 for the most part (see 
Sect. 14.2.2(1 . For f cS (M n ) we thus obtain 



(24) 



/eff = 7y\T 



— 1 dy{l - sin 2 a) = 1 - sin 2 a cff 
lY Jy, r 



where we used the midpoint rule. The angle a g can be 
interpreted as an average bending angle. As the ansatz for 
the Mach-number dependence of f e g we thus take 



eff 



1 — sin a c 



1 



V3 



M? 3 



(25) 



3.2.3. Energy dissipation 



A first expression for the column-integrated dissipated en- 
ergy per time can be obtained from energy conservation, 
£diss = £drv — £kin- For £drv we just derived an expression, 
Eas. 1211 and 1231 For £kin we obtain, within the frame of 
self-similarity, 



(26) 



• _ p m V 2 ms d£ cdi _ 3 V2 M 3- 01 

fckm- 2 dt -Pua 2 M U , 

where we used Eqs.dlEl and EE] to EH Together we get 

fdiss = PuaHll [1 - m MP 3 - 0.5 nlM-^]. (27) 

The energy dissipated per time within an average column 
of length £ c di is thus independent of this length. If en- 
ergy dissipation occurs only (as within the frame of Euler 
equations) or at least dominantly in shocks, which implies 
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that the average distance between shocks increases and 
/ or the average strength of the shocks decreases as the 
CDL grows. 

A second expression for £diss can be obtained from 
dimensional considerations. The energy dissipated per 
unit volume per unit time must be proportional to 
PdissWdiss^diss- Here, and ^diss are the character- 

istic density, velocity, and length scale of the dissipation. 
The energy dissipation within an average column of length 
4di can thus be written as £ d iss oc pdissVdiss^diss^dl- As all 
length scales must evolve equally with time within the 
frame of self-similarity, ^ c di/^diss must be constant, thus 



c 

£-diss OC /?diss^, 



iss^diss- 



(28) 



Comparison of Eqs. E3 and EHl suggests Wdiss oc aM u 
and a more complicated Mach-number dependence for 
Pdiss- As u rms is the only velocity scale we have, it seems 
natural to assume that i>diss oc u rms . It then follows 



that 



or ft = 1 (and ft = 0). We note 
that iGammie fc Ostriken l[l996|) even found v^iss — Vrms 



for a ID case. 



3.3. Summary of expected scaling relations 

If a self-similar solution exists, we expect the following 
dependencies: 



Pxn 


= tji p u M^ = ?7ip u , 


(29) 


-/Vf mis 


= m M^ = 77r 1/2 M u , 


(30) 




= 4di/T = 2r?f 1 aM u , 


(31) 


^drv 


= p u a 3 M^l- m M^), 


(32) 


^kin 


= p u a 3 Af u 3 0.5 r?!, 


(33) 


^diss 


= p u a 3 M*(l- m M^ -0.5 % 2 ). 


(34) 



Note the differences to the ID solution: Eq.lJHlpredicts the 
CDL mean density to be independent of M u and «2d oc 
aM u , in contrast to pid oc A'/ 2 and Kid oc a/M n . 

In deriving the above relations, we made four basic as- 
sumptions: a) we have simple Mach-number dependencies 
of p m , Wrms, and / c ff, Eqs.^^^J andl2~5l b) the CDL den- 
sity and velocity are uncorrelated; c) we have high Mach- 
numbcrs in the sense that r^M^f" 1 >> 1 or M r 2 ms >> 1; 
d) fdiss oc u rms . 

In Sect. 0] we are going to check the validity of these 
assumptions and confront Eqs. [2H1 to U2] with numerically 
obtained values. We expect good agreement as long as 
Afrms >> 1, thus dissipation in shocks likely dominates, 
and as long as £ c di << Y. The 'Euler character' of the 
solution should prevail under these conditions. We also 
determine those quantities that cannot be derived analyt- 
ically. These are, on the one hand, the coefficients r/i and 
773, as well as the exponent ft. On the other hand, there 
are quantities for which we have no analytical expression 
at all, like the wiggling of the confining shocks, the as- 
sociated distribution of the angle a, or the Mach-number 
dependence of the length of the confining shocks. 



4. Numerical results 

We now present our numerical results. After a brief phe- 
nomenological description of the solution in Sect. 14.11 
we give quantitative results for initial conditions 1=0 in 
Sect. 14.21 Results for initial conditions 1=1 and 1=2 are 
given in Sect. 14.31 and asymmetric settings are briefly ad- 
dressed in Sect. 14.41 Discretization and domain studies are 
the topic of Sect. 14.51 

4.1. Brief phenomenological description 

We begin with a brief qualitative description of the CDL. 
As an example, the density structure of run R22_1.2.2 is 
shown in Fig. [3] for three different times. 

A first characteristic is the local bending of the confin- 
ing shocks. The spatial scale of these wiggles increases lin- 
early with time, as the CDL accumulates more and more 
matter and gets more and more extended. The inclina- 
tion of the wiggles with respect to the direction of the 
upstream flows decreases with increasing upstream Mach- 
number (see Sect. I4.2~2*|l . Occasionally, we observe a su- 
perimposed 'bending mode' (e.g. bottom panel in Fig. |3J), 
which in appearance is somewh at similar to the bending 
modes of the NTSI described bv lVishmacl l)l994|) . 

A second characteristic is the patchy appearance of the 
CDL. The turbulent interior is organized in filaments and 
patches, regions within which a flow variable remains more 
or less constant. The spatial extension of these patches 
increases as well as the CDL accumulates more and more 
matter. The flow variables clearly mirror the supersonic 
character of the turbulence: the contrast between high- 
density filaments and extended patches in Fig. [21 easily 
reaches two orders of magnitude, the root-mean-square 
velocity is well above sound, and the mean density is sub- 
stantially reduced compared to the ID case. Shocks within 
the CDL are ubiquitous. 

4.2. Settings without CDL at t = 

For symmetric settings, and if there is no CDL at time 
t = 0, we expect to see the self-similar relations we derived 
in Sect. 13. 2l We express the time evolution of the solution 
we express in terms of 



£{N) = N/N 



Pn 



cdl 



PuY 



(35) 



This function monotonically increases at about the same 
rate as the mean extension of the CDL, since p m ~ ?7iPu 
(Eq.EH). In fact, p m » 30p u fSect. l4~2~H and thus £(N) = 
60 corresponds to ^ c di ~ 2Yq. For the symmetric case we 
consider in this paper, £{N) is proportional to the elapsed 
time. Using Eq. ^|to express N, we can write 



£(N) = N/N 



2rp u v u 



2v„ 



(36) 



PuY Y 

and £(N) — 60 then corresponds to a time r = 30Yq/v u . 
Or, if we use v u « 5w rms fSect. l4"2~T)l and Y ~ £ c di/2 for 
£{N) = 60, we obtain r « 2>£ c di/v vms . 
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Unless otherwise stated, averages and best fits in this 
section are always taken over the interval 10 < £(N) < 70 
and over all runs without CDL at time t = 0. The interval 
was chosen such that initialization effects have died away 
and that domain effects do not matter yet (Sect. B~5.1fl . 

We mention here already that the two most extreme 
simulations in terms of M u , R5JX2.4 and R87_0.2.4, of- 
ten differ somewhat from the other simulations. In the 
case of R5_0.2.4, we ascribe the deviation to the only sub- 
sonic turbulence and the correlation of density and veloc- 
ity (Af rms w 0.9 and corr(p, v) w -0.4, see Sect. 14.2.1)) . In 
the case of R87_0.2.4, the shocks become sometimes too 
strongly inclined with respect to the computational grid 
to be properly resolved by our numerical grid fSect. l4~2~2*|) . 

4.2.1. CDL mean quantities and correlations 

We first turn to the correlation of p and v and the CDL 
mean quantities p m and M rms , Eas.l29landl30l One of our 
basic assumptions in deriving these self-similar relations, 
namely point b) that the CDL density and velocity are 
uncorrelated, we find confirmed by our simulations. For 
nearly all symmetric simulations without initial CDL and 
for 10 < 1{N) < 70, we have 0.1 > corr(p, v) > -0.1. 
The only exceptions are the three low Mach-number runs 
Rll_0.2.4, Rll_0.2.2, and R5_0.2.4 with correlations of 
about -0.2, -0.2, and -0.4 respectively. The top panel of 
Fig. 0] shows the time evolution of corr(p, v) for five se- 
lected runs that differ only in their upwind Mach-number, 
5 < M u < 90. 

In the middle and bottom panel of the the same figure, 
Pm/ Pu and M rms /M u are shown as a function of l(N) for 
the same runs. Two things are apparent. First, the ratios 
take similar values for all five runs, indicating that indeed 
ft w and 02 ~ 1 for the exponents in Eqs. and l3"01 
Second, the ratios are not constant with i(N), indicating 
that the numerical solution is indeed only approximately 
self-similar. We come back to this point in Sect. [31 

To determine optimum exponents /3i, i = 1,2, we 
rewrite Eqs. |2H and |20| as equations for 771 and 772 and 
minimize the variance o" 2 (^i) ■ Considering all data points 
within 10 < £(N) < 70 of all runs without a CDL at 
t = 0, we find the smallest variances for f3\ = and for 
P2 = 1- The corresponding means are p(r]i) « 28 and 
p{rj2) ~ 0.21. Although clearly identifiable, the minima of 
g are relatively shallow. Changing j3\ or (3 2 by ±0.1, or 
excluding the very low Mach-number case R5_0.2.4 (for 
which Afrms ~ 0.9) changes a by only about 5%. By re- 
peating the analysis but allowing for a linear dependence 
of r\\ on £(N), we obtain the same optimum values for j3\ 
and /?2 but with considerably smaller variance. As i(N) 
increases from 10 to 70, 771 rises by about 25% (from 25 to 
31), while 772 decreases by about 15% (from 0.22 to 0.19). 

Part of our assumption a), namely the simple Mach- 
number dependencies of p m and M rms , thus seems justi- 
fied. With 772 = 0.21, assumption c), rfeM^ >> 1, is also 



fulfilled for most of our simulations. An exception is again 
run R5_0.2.4, for which r)\M^ ss 1. 

In summary, the simulation results, p m ~ 28/? u and 
M lms 1=3 0.21M U , essentially confirm the expected rela- 
tions, Eqs. HO and ED] % /2 ??2 = 1, predicted by Eg. I2TJ1 
is fulfilled to within 10% at any given time. The mean 
density is (nearly) independent of M u . As expected, the 
solution is only approximately self-similar, M lms decreases 
by about 15% as i(N) increases from 10 to 70. 

4.2.2. Confining shocks 

The turbulence within the CDL is driven by the upstream 
flows. The confining shocks of the CDL affect this driving 
in two ways. The less inclined the shocks are on aver- 
age with respect to the direction of the upstream flows 
(smaller angle a e g in Eq. I24fl . the more kinetic energy 
survives shock passage and is available for driving the tur- 
bulence. The smaller the spatial scale on which the angle 
a varies, the smaller the scale on which the energy in- 
put changes. In the following, we analyze how these shock 
properties depend on M u and on £ c &\. 

For this purpose, we specify the following basic quan- 
tities. The discrete x-position of the left and right shocks, 
s\ and s r , defined for each discrete y-position yj as the two 
cell boundaries where the Mach-number drops for the first 
time from its upwind value M u to 0.8M U . We determine 
the average extension of the CDL, £ c di, as 

1 - 

4di= 7 ^[s r (^)-si (37) 
J 3=1 

The length of the left and right shocks, ^ s h,i and £ s h,r, are 
computed as 

J 

4h,i = - + (Vj - vi-i?] 1 ' 2 , (38) 

where J is the number of cells in y-direction, and i — l,r. 
We define the angle a\^(yj) as the angle between the x- 
axis and the tangent to the shock (see Fig. Its numer- 
ical computation is described in Appendix lAl To obtain a 
number distribution, we sort the values a^ r (j/j) e [0, 7r/2] 
into 60 bins. Finally, to obtain a measure for the scale 
on which the shocks are wiggled, we look at the auto- 
correlation functions Lj r , 

F = <[ S - l {y 3 )-s;]-[s- l {y 3+ y cm )-s- l ]> ^ 

where is the variance of the shock position s;, and 
< . > denotes the mean over all discrete position yj. For 
each time, wc determine 2/corro 

such that Ti(y COI1 - ) = 0.5. 
Averaging y corro over both shocks gives a mean auto- 
correlation length ^corr, 

4orr = ^ [Z/corr (s\) + y CO rr (s r )] • (40) 

A larger auto-correlation length £ corr then indicates that 
the shocks are wiggled on a larger spatial scale, but it 
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does not give the scale of the wiggles in absolute units 
(see below). 

All four quantities, CDL extension, number distribu- 
tion of angle a, shock length, and correlation length, are 
shown in Fig. |SJ The first panel of Fig. [S] shows the es- 
sentially linear growth of the CDL with £(N). The growth 
rate, however, slowly decreases with increasing £(N). The 
slope of a linear fit in the range 40 < £(N) < 70 is 
roughly 10% flatter than the slope obtained in the range 
10 < £(N) < 40. This fits with the slight increase in 
p m , observable in the middle panel of Fig. 0] The sec- 
ond panel of Fig. [S] shows that the average shock length 
£ s h = 0.5(^ s h.i + ^sh.r) is fairly constant with respect to 
£(N) but increases with M u . Assuming a dependence of 
the form £ sh = VshYMP° h , the variance a 2 (rj B b) becomes 
minimal for /3 s h = 0.8. As can be seen, the two runs 
R5 _0.2.4 and R87_0.2.4 again behave somewhat differently. 
If we neglect these two runs, /3 s h remains unchanged but 
a is reduced by about 40%. The third panel of Fig. 03 
shows that larger upwind Mach-numbers lead to less in- 
clined shocks with respect to the direction of the upstream 
flows (lower values of a). Shown is the number distribu- 
tion of a, averaged over 10 < £(N) < 70. Individual runs 
show a slight shift towards higher values of a as £(N) 
increases. This shift is, however, small compared to the 
effect of M u . The fourth panel of Fig. shows the auto- 
correlation length £ COIT . It not only depends on M u but 
is also proportional to £ c di- The best fit is found to be 
£ corr w 0.74diAC° 6 - Th e fifth panel of Fig.0shows £ co „ 
scaled with this best fit. From these scaling properties of 
^corr, we take that higher values of M u lead to smaller 
scale wiggling of the shocks with respect to £ cc ji- 

The absolute value of £ COTI clearly depends on the 
choice of the threshold value in our definition, T(y corl ) = 
0.5. Figure HJ illustrates the variation of T\ as a func- 
tion of y corr at the example of run R43J3.2.4. The top 
panel of Fig. H3 shows that the initially present sinusoidal 
wiggling of the confining shocks does not get lost un- 
til about £(N) — 15, which is rather late compared to 
the other runs. Mode-like signatures again appear around 
£(N) <; 50. Our data give, however, no clear answer to 
how typical and persistent such signatures are. A basic 
problem is that their wave length soon becomes compara- 
ble (within a factor of 2 or so) to the domain size in the 
y-direction, which may affect the signatures. From the bot- 
tom panel of Fig. EI on the other hand, it can be taken that 
Fi essentially decreases linearly from 1 to about 0.2. The 
other simulations show a similar behavior. Consequently, 
the above scaling properties of £ COYI should also be ob- 
tained if smaller threshold values are used, down to about 

r(ycorr) = 0.2. 

Figs. 01 and also allow some insight into why runs 
R5.0.2.4 and R87JX2.4 sometimes fit not so well. The 
third panel of Fig. [3] shows that our spatial resolution 
is barely sufficient for run R87_0.2.4, the largest upwind 
Mach-number we have considered. The number distribu- 
tion here peaks at around a ~ 0.1. In terms of discrete 
positions this means that the shock position changes by 



about 15 cells in the x-direction as one moves from yj to 
Dj+i- Run R5J3.2.4, on the other hand, may deviate just 
because of its low Mach-number. The turbulence within 
its CDL is subsonic, M rms ss 0.9; and with rfeM 2 1.1 
and corr(p, v) sa —0.4 (Fig.0J top panel), it violates two of 
the basic assumptions made when deriving the self-similar 
scaling laws in Sect. 13.21 

In summary, as M u increases, the bounding shocks be- 
come less inclined with respect to the direction of the up- 
stream flows (smaller a), the fraction of upstream kinetic 
energy that survives the passage through the bounding 
shocks increases, and the bounding shocks themselves are 
wiggled on progressively smaller scales (smaller £ COII /£ c d\- 

4.2.3. Energy balance 

Energy input into the CDL occurs only at its confining 
interfaces. Energy dissipation, on the other hand, occurs 
throughout the CDL volume. Nevertheless, according to 
the analysis in Sect. 13.21 both £drv and fdiss should be 
independent of the CDL extension if dissipation is only 
due to shocks and if £ c di is small compared to Y. The 
average distance between shocks must then increase and 
/ or the average strength of the shocks must decrease as 
the CDL grows. 

To determine £ d rv we must compute the driving effi- 
ciency / ft = fdrvA^ekin.u- The corresponding integral in 
Eq. 1241 is evaluated numerically, and the resulting driving 
efficiency is shown in the top panel of Fig. [7\ As can be 
seen, larger Mach-numbers lead to more efficient driving, 
and a smaller part of the upstream kinetic energy is ther- 
malized already at the confining shocks. The driving effi- 
ciency /off increases by about a factor of four between runs 
R5_0.2.4 and R87_0.2.4. Also noteworthy is that the abso- 
lute value of the driving power S^rv differs by more than 4 
orders of magnitude between runs R5_0.2.4 and R87J3.2.4. 
The best fit for the assumed Mach-number dependence 
(minimization of (7(773) in Ea. I25JI yields (3^ = —0.7. The 
corresponding values of 773 = (1 — f c s)M^ 7 are shown in 
the bottom panel of Fig. [7| From the figure we take that 
the second part of our assumption a), the simple Mach- 
number dependence of / e ff , seems justified. The figure also 
shows that / e ff, and thus the driving power £drv, is not 
strictly independent of £ c &\ but decreases with increasing 
£{N). Repeating the best fit analysis but allowing for a 
linear dependence of 773 on £{N) again leads to (3% — —0.7, 
while 773 changes from 3.1 to 3.6 as £(N) goes from 10 to 
70. The average value of 773 is 3.3. Omission of the extreme 
runs R5_0.2.4 and R87_0.2.4 does not change the result. 

We determine the dissipated energy as fdiss = £drv — 
^kin (Sect. I3.2.2fl . where £^ n is the change per time of 
the kinetic energy within an average column of the CDL, 
and fkin is directly from our simulation data. Figure [S] 
shows the numerically obtained value £diss (top panel) 
and the theoretically expected value (Eq. l3"4l £^ BS (mid- 
dle panel), both in units of JF ekln! u = PuV^, as well as 
the ratio of the two (bottom panel). For better display, 
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the theoretical value, which must not depend on i(N), is 
shown as a (constant) function of i(N). For the constants 
in Eg. 1341 we used the numerically obtained average values, 
7/3 = 3.3, /? 3 = -0.7, and 772 = 0.21. We used 773 = 2.75 
only for R5JX2.4, in accordance with the bottom panel 
of Fig. The numerically obtained value was smoothed 
for better display using a running mean with window size 
A£(N) — ±1. The effect of the smoothing is illustrated in 
Fig. 03 with the example of run Rll_0.2.4. 

Looking at the data of £diss and Sdrv, three points 
may be stressed. First, £diss (Fig. |H1 top panel) mirrors 
<?drv = -^ckin.u/cff (Fig. top panel), and the values usu- 
ally differ by less than 10%. This is not surprising. It im- 
plies, however, that for larger upstream Mach-numbers, a 
larger fraction of the upstream kinetic energy is thermal- 
ized only within the volume of CDL and not already at 
its confining shocks. For M u ^ 20, the energy dissipated 
within the CDL exceeds 50% of the upstream kinetic en- 
ergy (Fig. [HI top panel). 

Second, the bottom panel of Fig. [S] shows that £^ ss 
and £diss agree to within 10% most of the time. Given 
the wide range covered (5 orders of magnitude in £<iissj a 
factor of 20 in M u , and an increase by a factor of 7 in 
£(N)), we conclude that the self-similar solution gives a 
good estimate. 

Third, from the same figure it can be seen that £diss 
generally decreases, except for run R5_0.2.4. Excluding 
R5 _0.2.4, a linear fit to £diss/£diss yields a decrease of 10% 
as t(N) increases from 10 to 70. A similar fit to fdrv/^drv 
with £^ v = p u vl(l - 3.3Af" ' 7 ) yields an even slightly 
larger decrease of 13%. The net dissipation, £diss/£drv, in 
fact increases by 3%. Thus, as the CDL size increases, the 
absolute dissipation within an average column decreases 
while the net dissipation increases. 

In summary, the predicted scaling laws, Eas.r^tolHH 
are - within the range of applicability - essentially con- 
firmed by the simulations. The fraction of upstream ki- 
netic energy dissipated only within the CDL, and not at 
the confining shocks, thus increases with M u . Best-fit anal- 
ysis for the numerical constants yields f c g = 1—3.3 M~°- 7 . 
Both £drv and £diss decrease slightly with increasing £ c( u. 
The net dissipation rate £diss/£drv increases, but only 
slightly (3% increase as £{N) goes from 10 to 70.) 



On dimensional grounds, we can define two length 
scales based on volume properties of the turbulence, 

N -i/2 g m 

4 kl „ = ■ km , (41) 

^diss 

Nv 3 

£ Vrma = ipsa, (42) 

^diss 

where £kin = 2V Jv P v<1 ^ s ^ ne average column integrated 
kinetic energy density. Here V is again the 2D volume of 
the CDL, introduced in Sect. 12.3.11 The two scales are 
equal up to a numerical constant if the density and ve- 
locity are uncorrelated, in which case we can replace the 
average over the product pv 2 by the product of the aver- 
ages of p and v 2 , £ kin = £ c diPmV 2 ms = Nv 2 ms . As this is the 
case in most of our simulations we look at only one of the 
above quantities in the following, ^ Ckin , shown in the top 
panel of Fig. ^| For better display, as £ Ckia inherits the 
large time variability of £dissj it is smoothed in the same 
way as fdiss i n the bottom panel of Fig. |SJ 

Assuming a relation of the form £ Ckin = a Ckin ^ corr , we 
obtain optimal fits (minimum of cr 2 (a 0kin )) for a Ckin « 
1.3. The fits become only slightly better if a weak linear 
dependence of a Gkin on £(N) is allowed (13% change as 
£{N) goes from 10 to 70). ^e kin /^corr is shown in the middle 
panel of Fig. Looking directly at the dependence of 
£ Ckin on £ c di and M u , we find £ ekia oc £ c d\M~ 6 . This is the 
same dependence we found for £ COII in Sect. 14.2.21 £ Ckia 
scaled with this best fit is shown in the bottom panel of 

Fig. ma 

With increasing upstream Mach-number the charac- 
teristic length scale £ Ckin thus decreases with respect to 
the CDL extension. This is consistent with our observa- 
tion that for the same £ c di the interior of the CDL shows 
finer structuring (patches, filaments) for higher values of 
M u . Figure ITT1 illustrates this observation with the exam- 
ple of runs Rll_0.2.4 and R33_0.2.4. Shown in the figure is 
div(v), as the flow patterns, especially shocks, are better 
visible in this quantity than in density. 

In summary, our simulations show that the inherent 
length scale of the turbulence is proportional to the auto- 
correlation length of the confining shocks, independent of 
M u and £ c d\- With increasing M u , both length scales de- 
crease relative to the CDL extension, £ Skia /t c dl oc M~ a6 . 
The appearance of the CDL, the size of its patches and 
filaments, behaves similarly. 



4.2.4. Length scales of the turbulence 

In Sect. 14.2.21 we looked at the scaling properties of 
the confining shocks and pointed out that shorter auto- 
correlation lengths icon imply smaller-scale wiggling, thus 
smaller scale changes of the kinetic energy entering the 
CDL. In the following, we show that the interface based 
quantity £ C on is proportional to the length scale derived 
from the volume properties of the turbulence. We take this 
as evidence of the tight coupling between volume and in- 
terface properties, between the turbulence and its driving. 



4.3. Settings with CDL at t = 

We performed additional runs to study the influence of 
an initially present CDL. Figure IT^l illustrates the results 
for some selected quantities. Shown are all the runs we 
performed with initial condition 1=0 (no CDL at t = 0), 
1=1 (moderate CDL at t = 0), and 1=2 (massive CDL at 
t = 0). 

Comparison of the 1=1 and 1=0 curves in Fig.ll2lshows 
that an initially present CDL of moderate column density 
(N = 14 7Vo) soon develops characteristics similar to those 
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found in simulations without initial CDL. A quasi-steady 
state is reached for 1{N) <; 40. The 1=1 and 1=0 curves 
then agree to within about a factor of two for volume 
quantities like p m and M rms (first two panels in Fig. I12|) . 
Agreement seems slightly better for interface related quan- 
tities. For (1 — / c ff) M - 7 , shown in the third panel of 
Fig. El the 1=1 and 1=0 curves lie more or less on top of 
each other. The same is true for £ Ckiri /£ c( i\ M®- 6 , shown in 
the bottom panel of Fig. IT21 

The situation is slightly different for runs with an ini- 
tially rather massive CDL (1=2, with initially N = 56 No). 
Also in these simulations the CDL gets more and more 
turbulent. For all the quantities shown in Fig. El the 1=2 
curves approach the 1=1 and 1=0 curves. However, it takes 
these runs much longer to saturate. Only for £(N) > 240 
the curves finally seem to saturate, at similar values as the 
1=0 and 1=1 curves. That saturation does indeed occur 
around that time is also supported by Fig. El As can be 
seen, the average angle distribution of the confining shocks 
for run R22_2.2.2 first shifts to higher and higher values as 
£(N) increases. It then stagnates for the last two averaging 
periods, 190 < £(N) < 250 and 250 < £(N) < 310. 

In summary, we conclude that our symmetric simula- 
tions all end up in a similar quasi-steady final state. An 
initially present CDL only delays the development. The 
incoming flows also manage to generate (and sustain) a 
similar level of turbulence also within an initially massive 
CDL. 

4.4. Asymmetric cases 

We also computed a few asymmetric cases, where the two 
upwind Mach- numbers are different, M\ ^ M r . For the 
same reason as given in Sect. we expect the solution 
to only depend on M\ and M r . These dependencies are 
more complicated than those assumed in Sect. |3] as we 
now have two different upwind Mach-numbers. The simple 
dependencies of Sect. should, however, be recovered in 
the limit Mi -> M r . 

The basic physical reason for the more complicated de- 
pendencies on the upwind Mach-numbers lies in the strong 
back coupling between the turbulence within the CDL and 
the driving of this turbulence by the upwind flows. Our 
asymmetric simulations demonstrate clearly (much more 
clearly than the symmetric simulations) that the turbu- 
lence crucially affects the driving: although M\ and M r 
are strongly different, the corresponding driving efficien- 
cies are about equal, / e ffj « / e ff,r- Thus the efficiency does 
not depend primarily on the upwind flow. In fact, Fig. 1141 
shows that for both symmetric and asymmetric runs / e ff 
(averaged now over both shocks) can be described well by 

/ off = 1 - M rm ° s 6 . (43) 

The angle distribution of the two shocks behaves accord- 
ingly in that it is similar for both shocks and determined 
by -Mrms rather than by either Mi or M r . 



A more detailed analysis of the asymmetric case, in- 
cluding an approximate analytical solution, will be pre- 
sented in a subsequent paper. 

4.5. Grid and domain studies 

The numerical results presented in Sect. l4.2l were all based 
on simulations with a domain Y = 4Yo and a discretiza- 
tion of 1.5 • 10~ 3 Y (R=2) or 2560 cells in the y-direction. 
Here we want to check whether these choices have any 
systematical effect on the numerical results of Sect. 14.21 

4.5.1. Different y-extent 

To check whether the size of the computational domain 
has any systematic effect on the results of Sect. 14.21 we 
performed some of the simulations again, but this time 
on smaller domains of Y = 2Yo and Y = Yo. We also 
performed one simulation on a larger domain Y = 6Yo. 

Figure 1151 illustrates our findings for simulations on 
domains Y = 2Yo and Y = 4Yq. M rms shows no system- 
atic effect and is, as such, representative of other volume- 
related quantities (Fig. El top panel) . As a typical rep- 
resentative for interface-related quantities, / e ff also shows 
no clear overall effect of the domain size (Fig. El middle 
panel). The quantity for which we find the most clear ef- 
fect is the auto-correlation length £ COIT (Fig. El bottom 
panel). However, even for € corr the effect sets in only for 
two of the four runs and only for £(N) 30, i.e. once the 
CDL extension reaches about half the size of the smaller 
domain. For the numerical results in Sect. 14.21 £ c a\ ~ Y/2 
corresponds to £(N) = 60. We conclude that the y-extent 
of the computational domain has no apparent systematic 
effect on these results up to ^(^V) <, 30 and probably even 
up to £(N) < 60. 

A systematic effect of the computational domain on 
the numerical solution does become apparent if the sim- 
ulations are carried on much longer. One pair of runs, 
R22_0.2.2 and R22_0.2.6, were carried on much longer, 
till £(N) w 200. For this pair of runs, Fig. El shows the 
evolution of M rms for each run, as well as their ratio, 
-M rmSi 3 y /M rmS! y. The run on the smaller domain appar- 
ently shows a faster decay in M rms after £(N) « 100. From 
Fig.El we take that the behavior of this one pair of runs is 
most likely the rule, and not the exception. The top panel 
of Fig.lTTIshows £ COII , scaled, for all the symmetric runs we 
have performed and whose domain has a y-extent < 2Yo. 
The bottom panel of Fig. 1171 gives the same quantity for 
all the runs with a domain extention > 4Yq. Comparison 
of the two figures shows that runs on a domain < 2Yo 
saturate around ^ CO rrM°- 6 w 1.6Y . For runs on a domain 
> 4Y , 4orr reaches much higher values. 

4.5.2. Different discretization 

The results presented in Sect. 14.21 were all based on sim- 
ulations with a discretization of 1.5 • 10~ 3 Y (R=2) or 
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2560 cells in the y-direction. To check the effect of the dis- 
cretization on our results, we repeated several simulations 
with coarser and/or finer discretization. These simulations 
indeed reveal a systematic effect of the discretization on 
the values of average quantities. Nevertheless, the general 
properties of the solution, its approximate self-similarity 
and Mach-number dependences, remain unaltered. Only 
the numerical constants T)i are affected. The changes are, 
however, small when compared to the differences between 
the ID and 2D solution (for example, p m = r/ip u in 2D, 
while p m = M^p u in ID). 

We find that finer discretization generally leads to re- 
duced turbulence. Using finer meshes we obtained larger 
mean densities and lower values of M Tms , as shown in 
Fig. E| The driving efficiency gets lower and the shocks 
become more inclined with respect to the direction of the 
upstream flows, and the angle distribution is shifted to 
lower values. The characteristic length scale £ Ckiri remains 
about constant if taken in units of £ c di- 

A possible explanation for the reduction of turbulence 
(smaller M rms ) on finer grids could be the dominance of 
shocks for the energy dissipation in the CDL. On a coarser 
grid, the network of shocks within the CDL is less dense. 
The divergence plots shown in Fig. H3I illustrate this effect. 
A closer analysis of this idea is, however, beyond the scope 
of the present paper. 

We stress that, so far, no convergence has been reached 
in our discretization studies. Looking at the comparison 
of the three runs R22.0.1.4, R22.0.2.4, and R22_0.4.4 in 
Fig^Jshows that each reduction of the cell size by a factor 
of two leads to a reduction of about 20% in M rms . This 
indicates that the resolution of 2560 cells in y-direction 
in our standard runs (R*_0.2.4) and of 5120 cells in the 
y-direction in the refined runs is still not sufficient. This 
should be kept in mind when interpreting these results or 
any results on shock bound turbulent structures in 2D, let 
alone 3D. 

Also, no clear picture emerges regarding the deviation 
of M rms from the constant value predicted by Eq. A 
linear fit to M rms for 10 < £{N) < 70 yields -12% for 
run R22_0.2.4 and -23% for the two times coarser run 
R22J3.1.4. For runs R43_0.*.4, the grid dependence is the 
other way round: R43J3.2.4 shows a decrease of -25%, the 
twice coarser run R43_0.1.4 decreases by only -15%. 

5. Discussion 

We want to address four points in this section. First, we 
sketch possible reasons for the slight difference between 
the numerical solution and the relations we derived in 
Sect. |3] Second, we look once more at the driving of the 
turbulence and, in particular, the back-coupling between 
interface and volume properties. Third, we briefly consider 
our results in an astrophysical context, in particular with 
regard to molecular clouds. Finally, based on preliminary 
numerical results, we sketch the effect of some additional 
physics. 



5.1. The numerical solution versus the analytical 
solution 

In Sect. 13.21 we suggested that a self-similar solution to 
our 2D model problem may still exist for the limiting case 
where the system approaches infinity. The relations de- 
rived in that section give a reasonable estimate for the 
numerical results of Sect.0] However, while M Tms is con- 
stant in Sect. l3~2l the numerical simulations show a grad- 
ual decrease in M rms already for small CDLs, t c &\ Y/2 
(15% decrease of M rms as £{N) increases from 10 to 70, 
Sect. I4.2fl . We have no firm explanation for this differ- 
ence. We sketch three possible effects in the following, but 
stress that the available data do not allow us to clearly 
distinguish between them. 

A first, obvious reason could be the finite y-extent of 
the computational domain, Y. It sets an upper limit on 
the total energy input into the CDL, thus on the amount 
of mass within the CDL that can be driven. Once the CDL 
has accumulated too much mass, the driving per unit mass 
weakens and the turbulence starts to weaken. The spatial 
growth of the CDL slows down while the average density 
increases. The following considerations on time scales may 
illustrate this point further. 

An upper limit to the time at which Y starts to affect 
the solution is given by the time t y at which £ c< \\ = Y. 
At later times structures may still grow in the x-direction 
(up to i c d\ at most) but cannot grow any more in the 
y-direction (where Y sets an upper limit). For the runs 
in Sect. |OJ 4di = Y corresponds to £(N) w 120 or 
t y = 12Yo/f r ms- A lower limit for the decay time scale 
of the turbulence may be obtained as follows. For the 
case of uniformly driven isot hermal hydrodyn amic turbu- 
lence in a 3D periodic box, iMac Low! l)l999(l has shown 
that the typical decay time once the driving is turned 
off, to, and the initial driving wave length, A^rv, are re- 
lated by to sa Adrv/^rms- Assuming that this result also 
holds for our slab, that Adrv = Y, and that driving is 
turned off completely it follows that to ~ Y/w rms , or 
to w 4Yo/v Tms for the runs in Sect. FOl However, driv- 
ing continues in our simulations and so the effective de- 
cay time scale of the turbulence is likely to be much 
longer than to. Finally, for the runs in Sect. 14.21 and 
a typical integration time of £{N) — 60 corresponds to 
about t = 6Yo/w rms , a typical turbulent crossing time at 
£{N) = 60) is T cross = 4di/v rms ~ 2Y /u rms . Comparing 
these different time scales makes it seem likely that at 
£{N) = 60, turbulence in the center of the CDL is still 
essentially driven, not essentially decaying. 

Our simulation data do not allow us to either clearly 
confirm or reject the hypothesis that the finite y-extent of 
the computational domain is responsible for the slight de- 
crease in Mrms that we observe at early times, £(N) 70. 
If the finite domain size were responsible, M rms should de- 
cay differently on different domains. Comparison of simu- 
lations on different domains up to £(N) w 70 fSect. l4~5~T|) . 
however, gives no clear picture. The data are rather noisy, 
and simulations on domains 2Yq and 4Yq show no system- 
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atic differences as long as £{N) <, 30 (^ c( ji < Y/2 on the 
smaller domain). Only for much later times, £{N) >> 70, 
well beyond the range for the results in Sect. 14.21 does Y 
have a clear effect and M rms decreases faster on smaller 
domains (Fig.lTfijl. 

A second, more speculative, reason might be numerical 
dissipation, provided that its effect were to increase with 
£ c di- While we have no evidence that the latter is really 
the case, it may also be hasty to disca rd this possibility 
right awav. IPorter fc Woodwardl ll 19941) found, by observ- 
ing how simple 2D hydrodynamical flows (shear flows and 
sound waves of definite wave number, their section 3.3) 
damp with time, that the decay rate due to numerical dis- 
sipation alone is a non-linear function of the wave number. 
Their results are certainly not directly applicable to the 
present case. But in view of these results, and given the 
change in structure size with £ c< u as suggested by Fig. [31 it 
might be possible that the effect of numerical dissipation 
indeed changes with £ c ai- Note that this would also imply 
that the MILES approach, outlined in Sect. 12.11 were not 
strictly valid for the problem we consider. The currently 
available data do not allow us to clearly reject or confirm 
the effect. 

A third reason, or rather an amplifying mechanism, 
could be back-coupling between M rms and the driving ef- 
ficiency. Once the turbulence within the CDL is slightly 
reduced (for whatever reason), the reduction is further 
amplified by the back-coupling between turbulence and 
driving, f e g = (1 — M~° a 6 ). The decrease in M rms results 
in larger inclination of the shocks with respect to the di- 
rection of the upstream flows, more energy is dissipated at 
the confining shocks of the CDL, and less driving energy 
enters the CDL. For the observed 15% reduction of M rms , 
the reduced driving may, in fact, play a dominant role: 
as £(N) increases from 10 to 70, £drv/£drv decreases by 
13% (Sect. l4~2.3(l . But to really estimate the relative im- 
portance of the three effects just sketched, further studies 
are certainly necessary. 

Two more points seem noteworthy to us in this section. 
One concerns the near independence of £diss on £ c d\- From 
Fig. El (increase in structure size with increasing £ c d\), we 
take that it is rather the increasing average distance be- 
tween shocks that allows fdiss to be essentially indepen- 
dent of £ c di and not so much the, on average, decreasing 
strength of shocks (Sect. 13.2.31) . Whether this is indeed 
true, only a closer analysis of the structure within the CDL 
along the lines of iMac Low &: OssenkopJ l)2000|) can tell, 
which is, however, beyond the scope of the present paper. 
Such an analysis could also shed light on whether (or in 
which sense) £ ekin (see Sect. I4.23|) is indeed a measure of 
the average distance between shocks. It would also allow 
us to quantify our impression that small scale structures 
are preferably located close to the co nfining interf a ces. I f 
true, this would fit with the result bv lSmith et al. 
that the high-frequency part of the shock spectrum is lost 
most efficiently. 

The other point concerns run R5_0.2.4. With 
corr(p, v) ~ —0.4 M rms w 0.9, it violates two of the 



basic assumptions we made in Sect. 13.21 Its mean den- 
sity is close to the isothermal value for strong shocks, 
p m w 22p u w 0.9p u M^. Both Sdiss an d £drv increase with 
£ c d\- With these characteristics, R5_0.2.4 may mark the 
transition from compressible supersonic turbulence, the 
topic of this paper, to compressible subsonic turbulence. 

5.2. CDL and confining shocks: a coupled system 

The turbulence within the CDL is 'naturally driven' in 
the sense that we control neither what fraction of the to- 
tal upstream kinetic energy, p u M^, really enters the CDL 
nor the spatial scale on which this energy input varies. 
Both are directly determined by the confining shocks in- 
stead and indirectly depend on the system as a whole. 
The driving efficiency at each confining shock scales with 
-Mrms, even for situations where M\ ^ M r (see Sect. 14. 4|) . 
The auto-correlation length of the confining shocks and 
the characteristic length scale of the turbulence within 
the CDL are proportional to each other, both scaling as 
£ c diM~ - 6 . We take these facts as evidence that the CDL 
as a whole, its interface and volume properties, forms a 
tightly coupled, quasi-stationary, and self-regulating sys- 
tem. Back coupling between post shock flow and shock is 
also d escribed in other contexts, for example bv iFoglizzol 
l|2002|) for the case of Bondi-Hoyle accretion. 

An aspect that remained elusive in Sect. 01 is the spa- 
tial scale on which the energy input varies, the energy 
injection scale. To really tackle this issue, it would be nec- 
essary to analyze the energy spectrum of the CDL. This 
task requires, however, some caution because of the highly 
irregular boundary of the CDL, and we postpone it for 
the moment. Nevertheless, we would like to present a few 
thoughts on the subject. 

A first question is whether it is justified to speak at 
all of only one injection scale, of monochromatic driv- 
ing. The homogeneous upstream flow is modulated by the 
confining shocks. These are wiggled on a variety of spa- 
tial scales at any given moment. This strongly suggests 
that the kinetic energy input into the CDL is most likely 
not monochromatic but occurs at a whole spectral range 
instead. Consequences of such non-m onochromatic driv- 
ing ha ve been studied, for example, bv lNorman fc Ferraral 
(Il99ffl . 

It also seems worthwhile to briefly look at 

monochromatically-driven tur bulence, in particular 

at the numerical simulations bv lMac Low! 1^999). For the 
case of artificially, monochromatically driven hydrody- 
namic turbulence in a 3D box with periodic boundaries, 
he found that the characteristic length of the turbulence 
is proportional to the driving wave length, independent of 
the Mach-number: A/^j™. — 1-42, where A is the (known) 
driving wave length and £^ d is the 3D analogon of £ c , . 
in Eq. E] In addition, IMac Low! ((1999) observed that 
increases with A, which is mirrored in the apparent 
increase in the structure size (patches, filaments). 
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Although our setting clearly differs from that 
of lMac Lowl l)l999|L two thoughts come to mind. The first 
is an actual observation, namely that we also observe an 
increase in structure size with . The second tho ught is 
more of a question or speculation. |Mac Low! (Il999h deter- 
mines the proportionality constant between the character- 
istic scale of the turbulence and the monochromatic driv- 
ing wave length. One may wonder about the implications 
of this finding if not one driving wave length is present 
but a whole spectrum. How will the characteristic length 
scale of the turbulence, which can still be determined fol- 
lowing Eq. 1411 depend on this spectrum? And, given our 



finding that £ e 



r , what does £ corr tell us about this 



spectrum? Both questions should become tractable once 
the energy spectrum of the CDL is determined. 

5.3. A glimpse at astrophysics 

With regard to astrophysics, the presented work basically 
suggests that, within the frame of isothermal hydrody- 
namics and a roughly plane parallel setting, larger Mach- 
numbers of the colliding flows results in a finer and finer 
network of higher and higher density contrast within the 
interaction zone. In different types of wind-driven struc- 
tures, this connection between Mach-number and struc- 
ture may be directly observable. 

For the clumping of line-driven hot-star winds, our re- 
sults suggest that the sheets or clumps formed by the in- 
stability of the line-driving are not homogeneous but pos- 
sess fine-scale substructure with high density contrast. 

Concerning molecular clouds, we first mention that re- 
cent ar guments suppo rt th e idea, origina lly brought for- 
ward bv lHunterl (|1979f ) and lLarson ( 1981J), that molecular 
cloud s result from th e colli sion of large-scale flows in the 
ISM. iBasu fc Muralil l)200ll) make the point that small- 
scale driving (« 0.1 ■ 1 pc) of molecular clouds is incom- 
patible with observed total luminosities, unless the energy 
dissipation rates derived from MHD simulations are seri- 
ously overestimated. Using a principal co mponent analysis 
of 12 CO (J=l-0) emission, iBruntl l|2003[) identifies large- 
scale flows of atomic material in which the globally tur- 
bulent molecular clouds are e mbedded. Similar observa- 
tional r esults were reported bv iBallesteros-Paredes et all 
l|l999ah . 

Driven supersonic turbulence as a structuring 
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length of the tu r bulenc e , and thus the largest structure 
size (jMac Low! Il999t IBallesteros-Paredes fc Mac Lowl 
2002), is usually a free parameter. Our results show 



instead that, at least for the case of an isothermal, 
shock compressed, supersonically turbulent 2D slab, the 
structure size rather depends on the size of the slab or 
cloud. 

5.4. Additional physics: an outlook 

The model presented in this paper covers only some 
very basic physics. To obtain results with a more 
direct relation to reality, additional physics must be 
included in the future, among these the following. 
Strongly asymmetric flows, where M\ ^ M r , lead to 
more complicated dependences, as we will demonstrate 
in a forthcoming paper. Inclusion of radiative cooling, 
instead of assuming isothermal conditions, can affect the 
problem in different ways. Therm al instability can lead 
to a d ditional dynamical e ffects llChevalier fc Imamural 



1982t iGaetz et all Il983 IStrickland fe BlondirJ I l995| 
Walder fc Fohnilll996l]Hennebelle fc Peraultlll99a bootj 



Vazquez-Semadeni et all l2~000t IKovama, fe Inutsukal 
2002 1 lAudit fc Hennebelle! 1200.4 iHeitsch et all 120051: 



Pittard et all 120051 iMiuiioiuJ l2l)()5li. Extended cool- 
ing layers, on the other h and, tend to a ct as a 
cush ion. Simulation s by IWalder fc Folinil l|l999ft 
and IWalder fc Fofinl fcOOObll which include radiative 
cooling but have otherwise similar parameters as some of 
the simulations presented here, show comparatively more 
small scale structure and even roll-ups at the interfaces 
confining the CDL. The CDL as a whole evolves less 
violently, and mean densities are about a factor of four to 
eight higher that what we found here for the isothermal 
case. Strongly asymmetric flows, where M\ ^ M r , lead 
to a qual itatively different solut ion if radiative cooling is 
included II Walder fc Folinilll 998ft and to more complicated 
dependences on the upwind Mach-numbers in the isother- 
mal case, as we will demonstrate in a forthcoming paper. 
The role of thermal conduction has only been considered 



by re l atively few publication s so far Bcgelman & McKcc 



Kovama fc Inutsuka 



I l990t iMvasnikov fc Zhekovl Il99i 

2004). Global bending of the interaction zone affects the 
stability properties of the interaction zone as a whole 
and thus probably also its interior properties. In colliding 
wind binaries, for example, matter is transported out 
of the central part of the system and diluted in the 
outer part. Simulations of bow shocks and colliding 
winds in binaries show strong traveling waves, together 
with a systematic change of the me an properties in the 
flow off from the stagnation point llStevens et all Il992t 
IWalder fc FoliTlll 995UBlondin fc Koerwerlll 998ft. 

6. Summary and conclusions 

We looked at symmetric, supersonic (5 ^ M u <, 90), 
isothermal, plane-parallel, colliding flows in 2D. The re- 
sulting shock-confined interaction zone (CDL) is super- 
sonically turbulent (1 <, M lma ^ 10). We investigated the 
CDL and its interplay with the upstream flows by dimen- 
sional analysis and numerical simulations. The latter we 
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generally stopped when £ c di ~ Y/2. The results are in- 
teresting not only with regard to flow collisions, but also 
shed new light on the properties of supersonic turbulence 
in general. 

The numerical simulations show that the CDL has an 
irregular shape and a patchy, supersonically turbulent in- 
terior. The driving of the turbulence is natural in that it 
depends on the shape of the confining shocks. The dimen- 
sional analysis is based on isothermal Euler equations in 
infinite space. Within this frame, a self-similar solution 
may exist that would depend on M n but must not depend 
on £ c dh Relations for average quantities are obtained un- 
der some further simplifying assumptions (SectEU- 

Based on both the analytical and numerical results, we 
arrive at the following conclusions. 

1) Comparison of the numerical and the self-similar 
solution shows generally good agreement if M rms 1. 
The modest deviation between the numerical and the self- 
similar solutions increases with l c( ±\. We suggest some ex- 
planations for the deviation, but our data do not allow any 
clear conclusions on the issue. For M rms <, 1, we have but 
one simulation. It shows clear differences to the other runs 
and may be more characteristic of compressible subsonic 
turbulence than of supersonic turbulence. 

— 1/2 

2) The CDL is characterized by M lma w ?] 1 M u and 
pm ~ ViPu- The average compression ratio of the CDL 
is thus independent of M u . This is in sharp contrast to 
the ID case, where p m! id = M^p u . From the numerical 
simulations, we find r/i s» 30. 

3) The turbulence within the CDL and the driving ef- 
ficiency are related by / c ff = 1 — Mj: 6 . The relation also 
holds for asymmetric settings, where M\ ^ M T , empha- 
sizing the mutual coupling between volume and interface 
properties. For larger upstream Mach-numbers, the shocks 
confining the interaction zone are less strongly inclined 
with respect to the direction of the upstream flows. The 
driving is more efficient, a larger fraction of the upstream 
kinetic energy is dissipated only within the CDL and not 
already at the confining shocks. 

4) The characteristic length scale of the turbulence, 
^e kln j and the auto-correlation length of the confining 
shocks, £ C orr, are proportional to each other. Both scale 
as £ c d\M~ Q - 6 , this although the former is based on vol- 
ume quantities while the latter is derived from interface 
properties. 

5) The separation of filaments and the size of patches 
within the CDL both get larger as 4di increases and/or 
M u decreases. 

For increasing upstream Mach-numbers in summary 
we thus expect a faster expanding CDL with less strongly 
inclined confining interfaces with respect to the direction 
of the upstream flows, similar mean density, finer interior 
structure relative to the CDL size, and a gradual shift of 
the energy dissipation from the confining shocks to inter- 
nal shocks within the CDL. We expect to observe these 
general dependencies in real objects where shock-confined 



slabs play a role, like molecular clouds, wind driven struc- 
tures, supernova remnants, or 7-ray bursts. 
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Appendix A: Numerical computation of 
obliqueness angle 

While shocks are smeared over approximately 3 grid cells 
in our simulations, the confining shocks in our analysis are 
specified as a series of discrete x,y-coordinate pairs only 
(see Sect. l4.2?2*|) . This information is sufficient to compute 
most shock-related quantities to good accuracy, for exam- 
ple the shock length £ s h. The only quantity that requires 
a more careful proceeding is the obliqueness angle a. If it 
were computed directly from the discrete shock positions, 
only discrete values would be obtained, for example 0°, 
45°, 63.4° etc. for one-sided differences. 

To compute the obliqueness angle ai(yj) (see Fig. ^ 
and Sect. I4.2.2|l at each position yj, 1 < j < J, of the 
left and right shock (si and s r ), we proceed as follows. 
In a first step, we use spline interpolation to double the 
number of points in the y-direction along the shock front. 
Next, we smooth the shock front slightly, using a running 
mean with an averaging window of ±5 points (this corre- 
sponds to an averaging window of ±2.5 points in the orig- 



inal data. Then we compute the derivative at each point 
of this smoothed shock front, using a 3-point Lagrangian 
interpolation. To avoid abrupt changes in the derivative 
from one point to the next, we smooth it again by a run- 
ning mean with averaging window ±5 points. We finally 
obtain the obliqueness angle ai(t/j), 1 < j < 2 J, as the 
arctan of the derivative. 

We checked that the size of the averaging window (±3 
points or ±7 points) has only a marginal effect on the 
angle distribution and the driving efficiency. For the latter, 
which is an integral over both shocks, tests show that a 
can even be computed directly from the discrete positions. 

Appendix B: List of runs, their parameters, and 
naming schemes 
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Fig. 3. The interaction zone of run R22_1.2.2, shown in 
density (logarithmic scale, in units of p u , color bar from 
to 4), for three different times: £{N) w 34 (top), £(N) w 
54 (middle), £(iV) « 74 (bottom). The spatial scale of 
patches, filaments, and wiggling of the confining shocks 
increases with £(N). 




Fig. 4. Time evolution of corr(p, v) (top), p m /pu (mid- 
dle), and M rms /M u (bottom) for runs R5.0.2.4 (dotted, 
dark blue), Rll_0.2.4 (dashed, purple), R22_0.2.4 (solid, 
red), R33_0.2.4 (dash-dotted, orange), R43_0.2.4 (dash- 
three-dots, green), and R87_0.2.4 (long dashes, pink). For 
these runs, £(N) = 60 corresponds to 4di ~ Y/2. 
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Fig. 5. Quantities related to the confining shocks: aver- 
age extension £ c( ji of the CDL (first panel), total normal- 
ized shock length I s d/(YM"' 8 ) (second panel), number dis- 
tribution (60 bins) of obliqueness angle a averaged over 
10 < i(N) < 70 (third panel), auto-correlation length 
^corr/Yo (fourth panel), and scaled auto-correlation length 
4orr/(4di-M~ a6 ), (fifth panel). Individual curves denote 
the same runs as in Fig. 0] 
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Fig. 6. Variation of T\, color coded, as a function of y cor i 
for run R43J3.2.4 (top panel). To allow for better display 
the color scale is limited to a range —0.5 < T\ < +0.5. 
Lower or higher values of T\ are uniformly colored in dark 
blue or red, respectively. For the same run, T\ is shown 
as a function of y corr for three selected times (bottom 
panel). £(N) = 30 (solid), £(N) = 50 (dotted), £(N) = 70 
(dashed). 
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Fig. 7. Driving efficiency (top panel) and best fit 773 
(1 — fctf)M®' 7 (bottom panel). For details see text. 
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Fig. 8. Numerically obtained (top panel) and theoreti- 
cally expected (middle panel) energy dissipation in units 
of the upstream kinetic energy flux density .F ekIlllU = PuV^. 
The constants in Eq. 02] were set to the best fit values, 
7/ 3 = 3.3, /3 3 = -0.7, and r] 2 = 0.21. We used ?y 3 = 2.75 
for run R5_0.2.4 (for details see text). The bottom panel 
shows the ratio of the two quantities. Individual curves 
denote the same runs as in Fig. ^ For better display, £diss 
was smoothed using a running mean with time window 
A£(N) = ±1. 
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Fig. 9. Effect of smoothing £diss with a running mean and 
window A£(N) = ±1, illustrated by run R33.0.2.4. Shown 



is Sdiss in units of T Ct 



before (dashed, black) 



and after (solid, red) smoothing, in units of erg cm s 
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Fig. 10. Characteristic length £ Ckin of the turbulence 
(top), in units of t cotT (middle), and scaled with best- 
fit e cd \M°- 6 (bottom) as functions of £(N). Individual 
curves denote the same runs as in Fig. 0] For better dis- 
play, £ c kin was smoothed by a running mean with window 
A£(N) = ±1. 
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Fig. 11. Plots of div(v) for two runs that are identi- 
cal except for their upstream Mach-number. Larger up- 
stream Mach-numbers lead, on average, to finer structure 
within the CDL and smaller scale wiggling of the confin- 
ing shocks. Shown are runs R33_0.2.4 (top) and Rll.0.2.4 
(bottom), both at a time when £ c( ji ~ 2Yo = Y/2. Blue 
(dark lines) indicates convergence, red (dark patches) di- 
vergence. 
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Fig. 12. Comparing runs with and without an initial 
CDL. Shown are p m /p u (first), M rms /M u (second), the 
scaled driving efficiency (1 — f c s)M®- 7 (third), and the 
scaled characteristic length of the turbulence £ 0kin /^cdi ■ 
M°- 6 for all symmetric runs. Line styles and colors denote 
initial conditions, (solid line, blue), 1 (dashed line, red), 
and 2 (dash-dotted line, orange). 
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Fig. 13. Time evolution of angle distribution for run 
R22_2.2.2. Shown is the average angle distribution for 
10 < £(N) < 70 (dashed, blue), 70 < £{N) < 130 
(dash-dotted, green), 130 < £(N) < 190 (dash-three- 
dots, orange), 190 < £(N) < 250 (long dashes, purple), 
250 < £(N) < 310 (solid, red). Also shown are the distri- 
butions for run R5_0.2.4 (black dots, right line) and for 
run Rll_0.2.4 (black dots, left line), both averaged over 
10 < £(N) < 70. 
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Fig. 14. Average / e g as a function of M rms for all our sym- 
metric simulations (triangles). In addition, we included 
data from our asymmetric runs (asterisks), for which 
1.6M r < Mi < 64M r and which initially have no CDL. 
Averages were taken over 10 < £(N) < 70 for simulations 
without initial CDL (blue triangles and green asterisks), 
over 40 < £(N) < 70 for runs with a moderate initial 
CDL (red triangles), and over 70 < £(N) < 140 for runs 
with a massive initial CDL (orange triangles) . Lines show 
f eS = 1-M| ms with ^ = -0.6 (dashed) and ^ = -0.6±0.1 
(dotted). 




Fig. 15. Comparing runs that differ only in the y- 
extent of the domain (Y = 2Yq and Y = 4Yo). 
Shown are M rmSj2 Y/M rmS! Y (top), / e ff,2Y//eflF,Y (mid- 
dle), and ^corr,2Y/^corrY (bottom). Individual curves de- 
note runs Rll_0.2.* (dashed, purple), R22.0.2.* (solid, 
red), R33.0.2.* (dash-dotted, orange), R43.0.2.* (dash- 
three-dots, green). 
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Fig. 16. Comparison of runs R22_0.2.2 and R22_0.2.6, 
illustrating the effect of a three-times larger y-extent 
of the computational domain on long time scales. 
Shown are M rms /M u (top) for R22_0.2.2 (solid, light 
blue) and R22_0.2.6 (dashed, dark red) and the ratio 
M rmS) 3Y/M rniS) Y (bottom). 
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Fig. 18. Comparison of M rms for runs whose spatial res- 
olution differs by a factor of 2 (subscript c = coarse, f 
= fine). Shown are (giving only the name of the finer 
run) runs R22_0.2.4 (solid, red), R22.0.4.4 (dash-three- 
dots, blue), R43_0.2.4 (long dashes, purple), and R11JX2.4 
(dash-dotted, orange). 
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Fig. 17. Scaled auto-correlation lengths of all symmetric 
simulations on domains with a y-extent less or equal to 
2Yq (top) and a y-extent greater or equal to 4Yo (bot- 
tom). 
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Fig. 19. Plots of div(-u) for two runs that are identical 
to run Rll_0.2.4, shown in Fig. 1111 except for their dis- 
cretization. The runs shown here were computed with two 
times lower (top) and at four times lower (bottom) reso- 
lution. Blue (dark lines) indicates convergence, red (dark 
patches) divergence. As can be seen, the number of con- 
vergent regions within an average CDL column decreases 
with decreasing resolution. 
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Table B.l. List of performed simulations. 



label 


M u 


t(N) 


^cdi/Y 




Pm 
P„ 


4h 

y 


/off 




Symmetric 


runs, no 


CDL at t=0 






R5_0.2.4 


0.42 




1.07 


n fin 

u.yu 


24 


1.1 


0.16 


Rll.0.2.4 


10.85 


88 


0.59 


2.2 


33 


1.5 


0.35 


R22.0.2.4 


21.7 


86 


0.30 


4.6 


30 


2.6 


0.59 


R33.0.2.4 


32.4 


86 


0.50 


6.9 


26 


3.6 


0.70 


R43.0.2.4 


A Q A 


QO 

OO 


0.55 


y.i 


on 

2y 


4.6 


0.76 


R87_0.2.4 


oO.o 


lUo 


0.82 


10. 


23 


12.1 


0.89 


R22.0.4.4 


21.7 


41 


0.25 


4.3 


35 


2.3 


0.55 


R22JU.4 


21.7 


88 


0.74 


5.0 


26 


2.7 


0.62 


R43JU.4 


43.4 


90 


0.59 


8.9 


32 


4.1 


0.76 


Rll.0.2.2 


21.7 


89 


1.10 


2.2 


33 


1.4 


0.33 


R22.0.2.2 


21.7 


307 


0.79 


4.7 


28 


2.6 


0.59 


R33.0.2.2 


32.4 


82 


1.45 


6.7 


30 


3.6 


0.70 


R43.0.2.2 


43.4 


73 


1.09 


9.4 


27 


4.7 


0.76 


R22.0.2.6 


21.7 


190 


0.84 


4.7 


29 


2.6 


0.60 


Symmetric runs, with CDL at t=0 


R22.1.2.2 


21.7 


87 


0.83 


3.3 


61 


1.9 


0.50 


R22.1.2.1 


21.7 


111 


1.33 


3.2 


68 


1.8 


0.46 


R22.1.4.4 


21.7 


199 


0.72 


3.4 


59 


1.9 


0.49 


R22.1.4.2 


21.7 


68 


0.40 


2.8 


91 


1.6 


0.39 


R22.1.1.2 


21.7 


115 


1.21 


3.9 


42 


2.4 


0.59 


R22.2.2.2 


21.7 


313 


1.44 


(2.4) 


(109) 


(1.5) 


(0.34) 


R22.2.4.2 


21.7 


186 


0.37 


(1.8) 


(253) 


(1.2) 


(0.24) 


R22.2.8.2 


21.7 


92 


0.14 


(1.4) 


(281) 


(1.2) 


(0.21) 



